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Abstract —  The  color  of  the  sea  is  determined  by  the  contents  of 
the  water,  especially  the  concentrations  of  suspended  particulate 
matter  (SPM),  phytoplankton  pigments  such  as  chlorophyll 
(CHL)  and  colored  dissolved  organic  matter  (CDOM).  Reversely, 
optical  sensors  that  measure  the  water-leaving  reflectance  spectra 
allow  us  to  calculate  the  desired  concentration  products.  In  this 
paper,  a  method  is  introduced  that  is  valid  for  both  case  1  and 
2  waters.  To  this  end,  model  is  fitted  to  reflectance  spectra, 
using  simulated  annealing  for  optimizing  the  mean  square  of 
the  reflectance  over  all  spectra. 

I.  Introduction 

A  variety  of  methods  exists  for  the  retrieval  of  water  quality 
variables.  A  first  is  based  on  the  decrease  of  the  blue  to  green 
reflectance  ratio  as  CHL  is  increased.  It  works  well  in  oceanic 
waters,  determined  by  phytoplankton  and  related  degradation 
products.  However,  it  fails  in  coastal  and  inland  waters,  where 
inorganic  suspended  matter  and  dissolved  and  particulate  non- 
phytoplanktonic  matter  also  affect  optical  properties.  Another 
method  calculates  CHL  using  the  red  to  near  infrared  ratio, 
which  is  effective  for  high  concentrations,  but  is  inaccurate 
for  low  concentrations.  The  multi-band  inversion  of  a  bio- 
optical  model  is  a  different  kind  of  method  that  has  found 
great  interest.  It  calculates  concentrations  for  SPM,  CHL 
and  CDOM  simultaneously  and  it  is  based  on  a  multi-band 
inversion  of  a  bio-optical  model.  The  forward  model,  such 
as  the  one  introduced  by  Kirk  [1]  and  Morel  and  Gentili  [2] 
relates  the  reflectance  just  below  the  water  surface  to  the  ab¬ 
sorption  and  back  scattering  coefficients.  Inversion  has  already 
been  proposed  in  the  literature  through  neural  networks  [3] 
or  matrix  inversion  [4],  The  matrix  inversion  is  restricted 
to  linear  forward  models  only,  in  which  case  an  analytical 
expression  can  be  found.  However,  these  linear  models  have 
constraints.  One  of  the  typical  problems  of  this  approach  is 
the  presence  of  negative  concentration  values  when  applying 
to  real  data.  Neural  networks  on  the  other  hand  have  to  be 
trained  with  large  data  sets.  Collecting  field  data  at  sea  is  time 
consuming  and  therefore  extremely  expensive.  As  a  result, 
existing  databases  are  rather  sparse.  An  alternative  apprach 
is  to  use  simulated  data  for  training,  at  the  risk  of  poor 
generalization. 

Here  concentration  values  are  estimated  by  fitting  a  model  to 


reflectance  spectra.  Some  optimization  schedule  and  criterion 
is  needed.  We  have  used  a  simulated  annealing  technique  (SA), 
optimizing  the  mean  square  of  the  reflectance  error  over  all 
spectra. 

II.  Analytical  model  for  the  irradiance 

REFLECTANCE 

The  subsurface  reflectance  can  be  related  to  the  remote 
sensing  reflectance  measured  from  (far)  above  the  water  sur¬ 
face  [5].  Through  models,  the  reflectance  is  related  to  the 
backscattering  and  absorption  coefficients,  from  which  water 
quality  variables  can  be  derived.  A  complete  derivation  of  the 
analytical  model  for  the  irradiance  reflectance  can  be  found 
in  [6]. 

For  simulations  of  subsurface  reflectance,  we  consider, 
throughout  this  paper,  water  bodies  which  consists  of  four 
optically  active  components:  pure  seawater  (w),  phytoplankton 
(ip),  non-algal  particles  (NAP)  and  CDOM.  Also  for  compar¬ 
ison  with  matrix  inversion  results  (see  section  IV),  we  use  a 
linearized  10P  model. 

The  total  absorption  coefficient  (a)  is  a  sum  of  contributions 
from  the  four  components: 

a(A)  =  aw(\)  +  a^(  A)  +  onap(A)  +  (Icdom(A)  (1) 

The  water  absorption  coefficient  aw  ( A)  was  taken  from 
Pope  and  Fry  [7],  Phytoplankton  absorption  coefficient  a^( A) 
is  computed  as: 

a0(A)  =BricA*CHL,  (2) 

with  Brie  A  taken  from  [8].  The  non-algal  particle  absorption 
onap  is  assumed  as: 

H'Nap(A)  =  unap(443)  exp(— Snap(A  —  443))  (3) 

with  Snap  =  0.0116nm_1  for  the  North  Sea  [9]  and 
Onap  (443)  assumed  to  be  proportional  to  the  suspended 
particulate  matter  concentration  (SPM)  with  the  constant  of 
proportionality  0.033.  Assuming  that  SPMnap  =  0.83SPM  [9], 
equation  (3)  can  be  written  as: 

0  033 

aNAp(A)  =  —  x  SPMNAp  x  exp(-0.0116(A  -  443))  (4) 
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The  CDOM  absorption  has  a  spectral  dependence  similar  to  for  scattering  coefficients.  We  assume  the  relationship  between 
Onap(A):  the  quantities  as  follows: 


ocdom(A)  —  acDOM(443)  exp(— S'cdom(A  —  443)),  (5) 

where  S'cdom  =  0.0167nm_1  for  the  North  Sea  [9]. 

Since  the  CDOM  scattering  is  negligable,  the  total  volume 
scattering  function  (VSF)  can  be  written  as: 

/3(A,  9)  =  /?u,(A,  8)  +  /3<jj(A,  9)  +  /3nap(A,  9),  (6) 

where  9  is  the  scattering  angle.  It  is  known  that  (3(f>{X19 ) 
is  highly  variable  depending  on  phytoplankton  species  and 
only  weakly  contributes  to  water  reflectance  (due  to  weak 
backscattering  ratio).  Therefore  instead  of  using  9)  the 

scattering  model  of  phytoplankton  and  covarying  (associated) 
particles,  (3c{X,9 )  was  used.  Defining  SPMnc  to  be  a  part 
of  SPM  which  is  not  covarying  with  CHL,  the  VSF  is 
decomposed  as: 

A,  9)  =  (3W{X1 9)  +  /3c(A,  9)  +  /3nc(A,  9),  (7) 


This  equation  can  be  rewritten  in  terms  of  scattering  coeffi¬ 
cients  6(A)  and  scattering  phase  functions  P(9): 

b{\)P{8)  =  bw(X)Pw(9)  +  bc(X)Pc(0)  +  6nc(A)Pnc(6)  (8) 


For  the  simulation  bw(X)  was  taken  from  Smith  and 
Baker  [10].  The  scattering  coefficient  due  to  phytoplankton 
and  covarying  particles  6C(A)  is  given  as  a  function  of  CHL: 

be  =  6C(660)  x 

with  6c(660)  =  0.407  x  CHL.  The  scattering  phase  function 
Pc(9)  is  computed  as  a  weighted  combination  of  two  Fournier- 
Forand  phase  functions  so  that  the  backscattering  ratio  should 
be 

=  6fec  (660)  =  0.0096  (10) 

be 

Details  on  how  to  make  this  scattering  phase  function  can  be 
found  in  [11]. 

The  scattering  coefficient  due  to  non-covarying  particles 
6nc(A)  is  modeled  as: 

000 

It  is  assumued  that  6nc  is  proportional  to  SPMnc  with 
^nc (A) /SPMnc  =  0.54  [12],  With  SPMC  the  dry  weight  of 
phytoplankton  and  covarying  particles  and  n  =  0.4,  equa¬ 
tion  (11)  can  be  rewritten  as: 

6Nc(A)  =  0.54  x  SPMnc  x 


The  scattering  phase  function  P^c(O)  is  computed  in  the  same 
way  as  for  Pc  (9)  but  using  a  constant  backscattering  ratio  of 
0.01833. 

According  to  the  IOP  model  described  above,  SPMnap  is 
needed  for  absorption  coefficients,  while  SPMnc  is  required 


SPM  =  SPMC  +  SPMnc  =  SPM,,  +  SPMNAp  (13) 
SPMC  =  0.234  x  CHL  (14) 
SPM0  =  0.5  x  SPMC  (15) 

Notice  that  SPMc  is  in  g/m3  and  CHL  is  in  mg/m3.  With 
the  equations  above,  the  three  unknown  variables  used  for  the 
simulations  are  CHL,  SPMNc  and  CDOM.  We  will  further  use 
the  term  SPM  for  SPMnc- 


III.  Methodology 


For  this  paper,  the  input  irradiance  reflectance  (R),  is  ob¬ 
tained  through  HYDROLIGHT  simulations  as  explained  in  IV- 
A.  For  real  imagery,  this  is  the  remote  sensing  reflectance 
measured  from  above  the  water  surface.  The  optimizer  will 
then  try  to  fit  a  modeled  reflectance  R  to  R.  The  variables  for 
the  optimizer  are  the  concentrations  CDOM,  CHL  and  SPM. 
We  used  the  model  from  (16).  Both  f  factors  are  implemented 
for  comparison.  The  first  depends  on  the  backscattering  ratio 
(with  bb(w)  the  backscatter  of  water)  [2]  and  thus  on  the 
concentrations  CDOM,  CHL  and  SPM.  The  latter  depends 
only  on  the  cosine  of  the  sun  zenith  angle  (/i o)  [13]  and  can 
therefore  directly  be  used  for  the  matrix  inversion. 
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/Kirk 
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0.63+  [  -0.22^—^ 

h 
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1+ 


/  bb{w) 

1  h 
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The  f  factor  can  be  recalculated  at  each  iteration,  using  the 
updated  values  needed  for  the  model.  Another  way  to  deal 
with  this  factor  is  to  treat  it  as  an  additional  parameter  to 
be  estimated  during  the  optimization.  This  option  has  some 
advantages  as  shown  in  IV. 

The  minimization  criterion  is  defined  in  matrix  notation  as: 

minimization  criterion  =  (R  -  R)  E"1  (r  —  R)* ,  (17) 

where  R  and  R  are  vectors  containing  the  measured  and 
modeled  reflectance  spectra  respectively  for  all  wavelengths 
and  S  is  the  covariance  matrix  of  the  signal  noise.  A  similar 
approach  was  followed  in  [14],  but  with  a  fixed  noise  value 
of  the  measured  reflectance.  We  have  adapted  this  formula  for 
a  band  selective  signal  to  noise  ratio.  When  this  is  known  for 
a  specific  sensor,  this  information  can  be  used  for  weighting 
the  criterion. 

In  general,  the  criterion  in  (17)  will  contain  many  local  min¬ 
ima.  The  global  minimum  will  thus  not  be  obtained  straight 
forward  from  an  analytical  result.  In  [14]  the  authors  used 
a  downhill  simplex  method.  In  this  paper,  we  used  simulated 
annealing  [15],  which  is  considered  more  reliable  for  finding 
global  minima. 


Simulated  annealing  was  also  applied  in  a  similar  context 
in  [16].  However,  the  authors  used  the  optimization  process 
only  once  to  obtain  a  unique  set  of  parameters  that  optimize  the 
semi-analytical  model.  In  this  paper,  the  optimization  is  used  to 
retrieve  the  concentration  values  for  each  individual  measured 
reflectance.  A  draw  back  of  this  method  is  that  it  is  time 
consuming  due  to  its  iterative  approach.  As  large  processing 
capacity  becomes  widely  available  through  low  priced  per- 
sonel  computers,  this  technique  is  no  longer  limited  to  small 
sized  test  samples,  but  can  well  be  used  for  real  imagery.  The 
processing  time  depends  on  the  spectral  resolution  (number  of 
bands)  of  the  signal.  On  a  Pentium  IV  single  2GHz  processor 
machine  an  image  containing  48  bands  was  processed  with  30 
pixels  per  second. 

Boundary  constraints  for  the  concentrations  are  easy  to  set, 
as  well  as  putting  initial  values  for  the  iterative  search  to 
the  optimum.  Constraining  and  intelligently  guessing  initial 
values  can  speed  up  the  process.  In  case  of  image  data,  the 
concentrations  of  the  previously  found  pixel  can  be  used  as 
initial  value  for  the  next,  as  nearby  pixels  will  likely  to  have 
similar  values. 

Using  hyperspectral  data,  the  optimization  process  for  fitting 
the  reflectance  spectra  is  constrained  by  a  large  number  of 
bands. 

The  proposed  methodology  is  not  limited  to  a  specific  sen¬ 
sor.  However,  hyperspectral  sensors  with  their  larger  number 
of  bands  will  constrain  the  optimization  process  more,  making 
it  more  robust.  Fewer  constraints  may  lead  to  ambiguous 
results. 

IV.  Experiments 
A.  HYDROLIGHT  simulations 

Experimental  results  are  obtained  using  a  monte  carlo 
simulation.  Random  concentrations  for  CDOM,  CHL  and  SPM 
are  used  as  input  for  HYDROLIGHT  4.2  [17],  a  state  of  the 
art  underwater  radiative  transfer  model.  The  obtained  dataset 
consists  of  subsurface  irradiance  reflectance  spectra  R(0— ) 
from  400  to  800nm  for  a  large  range  of  chlorophyll  (CHL), 
suspended  particulate  matter  (SPM)  and  colored  dissolved 
organic  matter  (CDOM)  concentrations  and  geometry  con¬ 
ditions  (sun  and  view  angles).  A  list  of  parameters  for  the 
HYDROLIGHT  simulation  is  summarized  in  Table  I.  The 
concentrations  are  randomly  variated  between  the  boundaries 
listed. 


TABLE  I 

HYDROLIGHT  PARAMETERS 


Parameter 

Value 

Solar  zenith  angle 

0°. 15°, 30°, 60° 

View  zenith  angle 

0°,  30.38°,  60.13° 

Relative  azimuth  angle 

0°, 90°, 180° 

Wind  speed 

5  m/s 

Sky  condition 

clear  (model  of  Harrison  and  Coombes) 

CDOM 

0-1  g/m3 

CHL 

0-30  mg/m3 

SPM 

0-100  1/m3 

Troughout  the  simulations  optically  deep  water,  which  is 
believed  true  for  turbid  waters,  and  homogeneous  water  col¬ 
umn  was  assumed.  In  the  current  dataset,  inelastic  scattering 
(such  as  chlorophyll  fluorescence  and  Raman  scattering)  and 
bottom  reflectance  was  not  considered. 

Published  IOP  models  based  on  measurements  in  the  North 
Sea  were  used  as  far  as  available.  Since  no  published  phyto¬ 
plankton  absorption  model  for  turbid  waters  exists,  an  oceanic 
water  model  for  phytoplankton  absorption  was  adopted.  De¬ 
tails  on  the  IOP  models  used  in  the  simulations  were  described 
earlier  in  section  II. 

B.  Results  and  Discussion 

The  concentrations  for  CDOM,  CHL  and  SPM  are  retrieved 
first  through  matrix  inversion  and  compared  to  the  simulating 
annealing  technique  (figures  1,2  and  3  respectively).  The  sun 
zenith  angle  is  used  for  calculating  the  f  value  according  to 
Kirk  in  (16).  Obviously,  the  Kirk  model  (represented  in  the 
figures  by  square  dots)  is  too  simplified  for  the  heterogenity 
of  angles  and  concentrations.  Moreover,  negative  values  are 
obtained,  ranging  from  -51  to  49  mg/m  '  for  CHL  concentra¬ 
tions.  For  SPM,  the  obtained  values  by  inversion  are  good  for 
low  concencentrations,  but  deteriorate  from  20  g/m3  onwards. 
Finally,  the  model  has  most  difficulties  with  concentrations  of 
CDOM. 

When  optimization  (SA)  is  applied,  the  f  value  can  be 
calculated  for  each  iteration,  as  explained  above.  For  this  study, 
it  has  been  calculated  using  Morel  and  Gentili  in  equation  (16). 
The  concentrations  obtained  for  CDOM,  CHL  and  SPM  are 
shown  in  the  figures  as  plus  signs.  The  results  are  much  more 
robust  for  different  sun  zenith  angles.  The  concentrations  are 
also  guaranteed  to  be  positive,  because  of  the  constraints  used. 

Finally,  figure  4  shows  the  modeled  irradiance  reflectance 
for  different  models  for  a  case  of  high  CHL  and  SPM.  The 
reflectance  spectrum  as  output  of  HYDROLIGHT  is  shown 
as  a  reference.  Simulated  annealing  has  been  applied  to 
these  models  to  fit  the  modeled  reflectance  irradiance  to  the 
simulated  reflectance.  The  concentrations  are  set  during  the 
optimization  process.  For  one  curve,  the  f  value  is  also  opti¬ 
mized  instead  of  calculating  via  one  of  the  equations  in  (16).  It 
is  shown  that  for  this  example  the  modeled  reflectance  is  closer 
to  the  true  (simulated)  reflectance.  Especially  the  retrieval  for 
SPM  will  improve  using  this  method.  This  is  part  of  ongoing 
research,  as  well  as  the  application  of  other  non-linear  models. 
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